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Abstract 

The singularities that arise in elliptic boundary value problems are treated locally by a 
singular function boundary integral method. This method extracts the leading singular 
coefficients from a series expansion that describes the local behavior of the 
singularity. The method is fitted into the framework of the widely used boundary 
element method (BEM), forming a hybrid technique, with the BEM computing the 
solution away from the singularity. Results of the hybrid technique are reported for 
the Motz problem and compared with the results of the standalone BEM and 
Galerkin/finite element method (GFEM). The comparison is made in terms of the 
total fiux (i.e. the capacitance in the case of electrostatic problems) on the Dirichlet 
boundary adjacent to the singularity, which is essentially the integral of the normal 
derivative of the solution. The hybrid method manages to reduce the error in the 
computed capacitance by a factor of 10, with respect to the BEM and GFEM. 
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1 Introduction 

Blamed for serious damages, in many engineering applications, singularities are under 
extensive computational investigation aiming to explore their origin and predicting 
their effects. A thorough review of several computational techniques that specialize in 
the treatment of the singularities of elliptic problems can be found in [1]. The cause of 
these singularities is found in, either a 'sudden' change in the boundary conditions (cf. 
the Motz problem in section 3), or the existence of a sharp comer in the geometry of 
the computational domain, both addressed in [1] and similarly treated. As for the 
abrupt changes in geometry, an example of practical interest would be the electric 
field singularities of conducting wedges surrounded by dielectrics and vice versa, in 
electrostatic problems [2]. A case of particular interest is the field singularity that 
arises in electro-capillary systems at the contact line (abrupt change in geometry), 
which is investigated in connection with phenomena that limit the electrostatically 
assisted wetting, such as the contact angle saturation [3]. 

A large family of techniques treating elliptic problems with singularities accounts for 
the asymptotic expansion of the singularity. In the present work we restrict ourselves 
to problems governed by the Laplace equation, which is the most common 
representative of the family of elliptic boundary value problems. The analysis of the 
singularity of the Laplace equation posed in a 2-d arbitrary domain - which will be 
the basis of the proposed method - produces an asymptotic expansion of the solution. 



u, that reads, u = 2^a,r^''f,(0)whQrQ, ai are the unknown singular coefficients (in 

fracture mechanics known as generalized stress intensity factors), r is the radial 
distance from the center of the singularity, 6 is the angle with reference to a boundary 
of the 2-d domain, jii and // are predetermined by the analysis of the singularity (cf 
section 2). 

The goal of this work is to embed the asymptotic expansion of the solution into a 
widely used computational method for analyzing physical systems governed by the 
Laplace equation. The method of choice is the boundary element method (BEM) [4, 
5] that is already a very commonly used technique in elasticity and potential problems 
and it is gaining ground in other types of problems, mainly due to its reduced 
computational cost, compared with other methods, e.g. the finite element method. 
Considering that matter, the greatest merit of BEM is that it reduces the dimension of 
the computational problem by one, i.e. only the boundary of the computational 
domain is discretized. 

The BEM has already been enhanced by singularity techniques; for example in [6-8] 
the BEM is augmented with singular boundary elements, i.e. elements with special 
basis functions that account for the singularity instead of the common polynomial 
basis fianctions. However, to embed as many leading terms of the asymptotic 
expansion of the singularity as the basis functions, requires the construction of many- 
node elements which in turn demands much tedious work. 

In this work the BEM is augmented with elements of the technique in [9], referred to 
as singular function boundary integral method (SFBIM), which employs the 
asymptotic expansion of u in the solution procedure, with as many singular terms as 
required. A similar treatment of singularities with a hybrid BEM was presented in 
[10]. The SFBIM was introduced in [9] and thereafter was efficiently applied to 
problems governed by the Laplace equation, such as the L-shaped domain problem 
[11]. Moreover, SFBIM was applied to problems governed by the biharmonic 
equation, such as the Newtonian stick-slip flow problem [12] and fracture problems 
with crack singularities [13]. 

The proposed method is essentially a coupling of BEM and SFBIM that results in a 
novel hybrid method and will be referred to as hybrid BEM/SFBIM. The 
effectiveness of the coupling lies in the similarities of the two methods, with both, 
being boundary integral methods that reduce the dimension of the computational 
problem. Briefiy, the coupling of the two methods is done as follows: The 
computational domain is decomposed in two subdomains, the first being a small 
segment of the original domain, surrounding the singular point, where the SFBIM is 
applied, and the second its complement, where the BEM is applied - the coupling is 
analyzed in detail in the following sections. 

2 The hybrid BEM /SFBIM 

Consider the Laplace equation posed in a 2-d arbitrary domain, Q., depicted in Fig. 1 . 
Let To and Fa? be the boundaries on which the Dirichlet and Neumann boundary 
conditions are imposed, respectively. The boundary of Q. is smooth except at the 
points O2, O3, where the boundary forms the angles ©2, ©3, respectively; Oi lies on a 



straight boundary segment with 0i=7i, therefore the boundary is smooth at Oi. 
Singularities arise at the points of the boundary at which a sharp corner is formed (O2) 
or there is a sudden change of the boundary condition (Oi) or both (O3). The change 
of the boundary condition can be between Dirichlet and Neumann or of the same 
kind, e.g. homogeneous-inhomogeneous Dirichlet. The asymptotic solution near Oi, 
O2 and O3 can be derived through separation of the independent variables, r, 6, where 
r is the radial distance from Oi and 6 is the angle with reference to the boundary. The 
boundary segments adjacent to 0\, O2 and O3 are straight lines in order to provide a 
simple form of the asymptotic solution. A more general expression for curved 
boundaries can be found in [1]. In this work, however, we restrict ourselves to the 
analysis of the simplest cases, that is homogeneous Dirichlet and/or Neumann 
conditions imposed on straight boundaries adjacent to the singular points. The general 
expression of the asymptotic solution for Neumann-Dirichlet (with ^ = on the 
Neumann boundary) conditions is 
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Fig. 1 . A 2-d arbitrary domain with Dirichlet (To), Neumann (Fm) boundary conditions and singularity 

source points Oi, O2 and O3. 



u ^UQ{r,0) + '^air^' cos(//,^) 



1=1 



and for Dirichlet-Dirichlet conditions 



u = Mo(r,^) + ^ sin(//,^) 



1=1 



where «/ are the unknown singular coefficients, jui are the known powers of the 
singularity that depend on 0, (0i, 02, or 03) and the types of boundary conditions (the 
details for the extraction of /ui are given in [1]), uq are particular solutions, which 
vanish for homogeneous boundary conditions, thus giving 



(1) 



1=1 



and 



(2) 



1=1 



with (1) being valid as well for homogeneous Neumann-Neumann conditions. 



As shown in Fig. 1, multiple singularities can exist in the same domain and we assert 
ourselves that the local expression of each singularity is valid as it was extracted 
neglecting the presence of the other singularities. The effect of a pair of singularities 
on the solution is investigated in section 3. 

The original domain of Fig. 1 is decomposed into three non-overlapping subdomains; 
the small subdomains Q.\, Q.2 contain the singularities - for simplicity we neglect the 
singularity at the point O3 ~ and the large subdomain Q.3, contains the bulk space of H 
where the effect of the singularity is relatively weak (Fig. 2). The subdomains are 
separated through a circular segment that surrounds the singular points at a given 
radius, R. The choice of the shape of the artificial inner boundary, that is being 
circular, is not justified by any other means than simplicity of implementation and the 
fact that is the most straightforward approach since the strength of the singularity 
depends on the radial distance fi-om Oi. The boundaries of each subdomain are, for D.\ 
r = , for Q2 r = and for Q3 T = Tj u r4 u u u - the small straight 

segments of the boundaries of Q.2 and Q.3 are not taken into account for the present 
analysis due to the properties of the suggested method, as it will be discussed below. 
The boundaries are grouped with reference to the type of the boundary condition as, 
= Fj , r^, = u Fj and the internal boundaries, F4, F5, don't fall into either one of 

the above types of boundaries as they do not have a specified boundary condition. The 
solution is approximated via (1) and (2) only in and ^2, respectively, while in Q3 
the solution is approximated via standard polynomial interpolation functions that are 
typically constant or linear. 



Fig. 2. Domain decomposition for the hybrid boundary integral method with two singularities; Fi = F^r, 



In we apply the standard BEM. In more detail, starting from the fundamental 
solution of the Laplace equation, the BEM extracts the following boundary integral 
equation 





n)u{^, n) + Jm(x, y) 



T 



dG{x,y;^,n) 
dn 



dr = j 



T 



du(x,y) 
dn 



G{x,y;^,n)dr 



(3) 



where G{x,y;^,n) is the free-space Green's function, defined as 
G(x,y;^,n) = ^ln(l/-J(x-^)^ +(>'~")^) ^nd n is in the direction of the unit vector 



normal to the boundary. The parameter X depends on (if, n); if {£,, n) e Q3 then A=l, if 
^ Q3 or 5Q3 then 1=0; if {B,,n) ^dO..^ cf below. The boundaries of Q3 are 

tessellated into a finite number of constant or linear elements, with one node placed in 
the center of each element or two nodes placed at the endpoints of each element, 
respectively. The solution, u(x,y), and its derivative, du(x,y)/dn, on the boundary, 

r, are approximated in terms of the basis functions, (x, y) , as 



uix,y) = J^u^jix,y), ^^^^^ = ^qj^' ix,y) (4) 



where Uj and qj are the nodal values of u and its derivative, respectively and is the 
total number of nodes. The basis functions, ^^{x,y), are piecewise polynomial 
functions that can either be constant or linear. In this work we use linear basis 
functions, unless indicated otherwise. 

By collocating the points (^, n) with the nodal positions (x/, yi) and inserting (4) in (3), 
a discretized version of (4) is derived that reads 



]=\ T ^ J=i r 



(5) 



The set of equations (5) can be written in matrix form 

Hu=Gq (6) 



where 



f on 

= A{x, , j;, ) + h ^ (X, y) f ^' ' ^- ^ dT, for i = j (7) 
f on 

Gy = ^0\x,y)G{x,y;x^,yi)dr 



The parameter A in (7) is equal to Vi if the z-th node lies in a smooth segment of the 
boundary, e.g. the middle of an element (that is the case for constant basis functions). 
When the i-ih node lies in a corner of the boundary, e.g. the endpoints of the element 
(that is the case for linear basis functions), 1 depends on the angle formed by the 
elements that are adjacent to the z-th node. It is convenient to circumvent the explicit 
computation of A with a simple technique without any loss of generality of the method 



[4]. This is done by applying a uniform u along the boundary, which in turn bounds 
the normal derivative of the solution, q, to be zero. Thus, (6) becomes 

Hu=0 

which in turn provides the diagonal elements of H using the computed off-diagonal 
elements, 

N 

instead of using the second equation of (7) that requires the computation of X. 

On each boundary of the domain there is a specified boundary condition, either 
Dirichlet where the uj are defined or Neumann where the qj are defined, excluding the 
internal boundaries where no boundary condition is applied. Thus, we can reorder (6) 
so that the unknowns he in the left hand side of the equation 

Ax = b (8) 

Based on the setup of the problem of Fig. 2, the matrix and the vectors of the 
unknowns are, respectively 

A = [ [H^ -[G],,,^ -[G],,,3 [h],,,^ [h],,,^ -[g] 
X = [ur, , qr, , qr, . Up , u^^ , q^ , qr^ f 

A is a Ny.(N + M^^ +^t, ) matrix, where M^^ and are the number of nodes on 

r4 and Ts, respectively. The extra Mp^ +^r, columns correspond to the two 

rightmost set of columns of A, - [g]^.^^ , - [g]^.^p_ , which are gathered in the LHS of 

(8) since r4 and Ts are internal boundaries and thus, both u and q are unknown. The 
extra + M^^ equations needed, will be provided by the coupling with the SFBIM 

that is applied on the subdomains, Q.i and Q.2- 

In Qi and Qithe solution, u, is approximated by (1) or (2) that can be rewritten 

u = j;^a,W'{r,0) (10) 

1=1 

where are harmonic functions of r and 6. The SFBIM incorporates the Na leading 
terms of (10), however, the proposed method performs sufficiently well with just a 
few leading terms as it will be seen in section 3. Thus, (10) is rewritten 



(9) 



1=1 



(11) 



The derivative of m normal to the boundaries of Qi and Q2 is approximated by 



t^ = £<.,m (12) 

on ,=1 on 

The Laplace equation is weighted with W'^ and Green's theorem is applied twice. The 
double integral that contains the term, V^W'' , vanishes since W'^ are harmonic and 
thus, the following boundary integral equation is obtained 

\—w'dr-\u^—dr = o k = \,...,N^ (13) 

Next, we apply (13) on the subdomains Q.\ and Q.2. The solution, u, is approximated 
on boundaries of the subdomains via (11). Its derivative, however, is approximated on 
the straight segments of fii and Q.2 via (12) and on the inner artificial boundaries, r4 
and T5, via polynomial basis functions (cf. Eq. (4)). This approach provides the 
equality constraints for the derivative of the solution at the artificial inner boundaries 
(C^ continuity constraints). The resulting set of integral equations from (13) reads 



7-k 

M^'rifi, 'rin, 5« ' { dn (14) 



The second boundary integral of (14) vanishes because either W or dWjdn is zero - 
depending on the boundary condition (homogeneous Neumann, dW/dn\ =0 and 

homogeneous Dirichlet, W\ = ) — on the straight boundary segments of Qi and D.2', 

the same applies for the third boundary integral of (14) for FctSQj. Thus, (14) 
becomes 

Y^qj lo^W'dr-Y,a, I W'^dT^O k = l,...,N 



7=1 Teen, '=1 Teen, 



dn " (15) 



Here are introduced along with the A^^ equations of (15), also Na unknowns (the 
leading singular coefficients) for each subdomain that contains a singularity; for the 
subdomains Qi and D.2 are introducedJV^ and A^^ j^^ unknowns, respectively. The 

rest of the constraints are provided by the matching requirement, that is to equalize 

the approximations of u, weighted by the basis functions O^, of the BEM and 
SFBIM along the boundaries, r4 and Fs (C° continuity constraints) 



^M. Jo^O't/r-^a, jw'O'dr = i = l,...,M (16) 

where M = Mp ,Mp^ and Mp ,Mp are the number of elements on r4 and Ts, 
respectively. From (16), M^^ + constraints are introduced to the problem, and 
overall, from (5), (15) and (16) we gather N + M^^ +^«,n, + ^«,n, equations 

with the same number of unknowns. The system (8) is completed with (15) and (16) 
and A, x are in expanded form (cf. (9) for the corresponding incomplete system) 



lH],.r, -[G],,r, -[G],,r, [h],,,^ [h],,,, -[g]„,^ - [g] ,,,, 

dn 



foWi/r - {W — dT 

J J Pin 



Jowi/r -jw—dr 

r, r, 



dn 

|a>3»cr -|wa>(/r o 
j<D(Ddr -jw<Ddr 



X = [ur, 'Qr, 'qr3 »r, »r, 'Qr, 'Qr, 'Wq, ^^n, } 



where a„ and a^^, are the vectors of the leading singular coefficients of the 
singularities that are contained in and ^^2, respectively, withaj-,^ =[«!,..., ] 
and ttf^^ =[«!,... ]. 



Numerical experiments 



The proposed technique is applied on the standard benchmark Motz problem [14, 15] 
(Fig. 3). The problem is governed by the Laplace equation posed in a rectangular 
domain, Q. = [-1,1] x [0,1] that is divided into five boundaries, denoted Fi,..., F5. The 
singularity is centered at the origin of the axes at the intersection of the boundaries Fi, 
F5, where there is a sudden change in the boundary condition from u\ = to 

du/dn\ = . On F3, F4, is applied the homogeneous Neumann boundary condition. 



du I dm 



and on F2 a Dirichlet boundary condition, u\ 



500. 




Fig. 3. The Motz problem domain; O = [-1,1] x [0,1]. 



The asymptotic solution of the singularity for the Motz problem is given by the 
infinite series 



cos 



1=1 



21-1 



e 



(17) 



which is exact for the entire domain and thus, true for any subdomain that includes O. 

Therefore, Q is decomposed into two non-overlapping subdomains as shown in Fig. 4 
where fli contains the singularity and contains the bulk space of the original 
domain. The harmonic functions for the Motz problem for the subdomain fli are 



W=r 



cos 



lk-\ 



e 



(18) 



which are valid for the entire domain; however, that is not a requirement for the 
proposed method. Only the leading A^^ harmonic functions are incorporated in the 
solution procedure, i.e. A: = 1,..., A'^^, where iV« is the number of the leading singular 
coefficients. 




r, r, o r, r, 

Fig. 4. Domain decomposition of the Motz problem; Q. = Q.\\J Q.2=\-\,\'\ x [0,1]. 

The application of the BEM on the domain Q2 depicted in Fig. 4 produces a system of 
equations as in (8) for a problem with arbitrary domain. The LHS reads 



(19) 



The system is augmented with N^+M equations resulting from the application of 
the continuity constraints on Q.\ of Fig. 4 



N 



N„ 



Y,Uj jo^O'cT -J^aijw'O'dT = i = l,...,M 



(20) 
(21) 



The resulting system from (19), (20) and (21) has a size N + M + N ^ vector of 
unknowns that reads x = [u^ , q^^ , u^^ , , q^^ , u^^ , qp^ , a]^ . 



The efficiency of the proposed method is evaluated in terms of the computed singular 
coefficients, «/, compared with their exact values found in [9]. Table 1 shows a set of 
results for the leading singular coefficients with typical solution parameters such as 
the discretization of the boundaries and the radius of Qi. The internal artificial 
boundary and the external boundaries of the Motz problem are discretized uniformly, 
e.g. if A^= 100 then each side of the rectangular domain has 50 equally sized elements 
and if M = 10 then the boundary Eg is represented by an even polygon with 10 sides; 
if linear or constant elements are incorporated the boundaries are represented as 
straight segments even though they may be curved. The results of the proposed 
method are in good agreement with the exact values of the two leading singular 
coefficients beyond which the discrepancies become so large that can possibly exceed 
several orders of magnitude - it should be noted that large deviations of the computed 
singular coefficients, fi-om the third and above, do not affect the computed uj or qj. In 
addition, regarding the solution in the subdomain where the SFBIM is applied, there 
is no 'noticeable' difference of the solution (that has the form of (11)) when the third 
leading singular coefficient and above is largely miscalculated - this is restricted to 
relatively small values of R with respect to the size of the computational domain, for 
this case -0.1. To further illustrate this point, we can apply an arbitrary value of as, 
a4, etc. and the deviation of computed solution (from (11)) from the exact (again from 
(11) using the exact «/) will not exceed margins of practical interest; this is the case 
when R is relatively small. When R is further increased, Na has to be increased as 
well, in order to achieve the same levels of accuracy (cf. below for a quantitative 
investigation (Fig. 5)). In problems with practical interest, a sensible choice of R 
requires only Na=\ or 2, since its value would be relatively small. However, to 
evaluate the performance of the proposed method Na is increased up to 10, while 
maintaining a small value ofR. Further increase of Na would be unnecessary since it 
has a little effect on the computed ai as seen in Table 2. This increase is limited 
however to a relatively small range of values of Na because otherwise the matrix A 
becomes ill-conditioned and the accuracy of the solution is compromised. 

Table 1 



Five leading singular coefficients of the Motz problem derived with the proposed method and their 
exact values; A'^ = 500 (total number of elements), M = 100 (number of elements on inner artificial 
boimdary), R =0.1, Na= 5. 







"2 


"3 


a4 


as 


Values of the 
proposed 
method 


401.067 


84.7409 


32.5103 


-153.517 


994.894 


Exact value 


401.162 


87.6559 


17.2379 


-8.07121 


1.44027 



Table 2 



Dependence of the a/ on A/;,; A/^ = 1 00, M = 20, = 0. 1 . 







Na=2 






Exact value 


«1 


401.240 


401.223 


401.231 


401.221 


401.162 


ai 




82.9425 


82.9647 


83.0262 


87.6559 


ai 






39.9397 


39.4761 


17.2379 


an 






-213.750 


-210.720 


-8.07121 


as 






1940.94 


1915.81 


1.44027 


«(, 








-11668.5 


0.331054 


«7 








158037 


0.275437 



«8 








-806343 


-0.0869329 










13477300 


0.0336048 


«10 








-65620700 


0.0153843 



The findings above suggest that a proper computational practice regarding Na is to 
incorporate only the leading singular coefficients that are needed (see below). 
Exceeding a reasonably small value of Na does not improve the solution, although it 
causes larger computational cost by further ill-conditioning the matrix A. 

The required value of Na depends on the quantity we wish to compute with high 
accuracy and also the radius, R. A very useful quantity for many engineering 
applications is the total flux on Dirichlet boundaries. For example, in electrostatics, 
the capacitance, which is traditionally used instead of the term total flux, is a quantity 
of primary interest [16, 17]. In this work and throughout the text we adopt the term 
capacitance (instead of total flux), denoted C, computed on boundaries with Dirichlet 
boundary conditions, e.g. Fg of Fig. 4. It is defined as 



J dn 



(22) 



If the boundary belongs to a subdomain that contains a singularity and specitically for 
the subdomain Q.i of the Motz problem of Fig. 4, then dujdn is approximated by (12) 
and (18). Thus, (22) on gives 



(23) 



With 



K, - sin 



2/-1 



-n 



^ 2M 
R 2 



(24) 



where Ki can be viewed as a weighting term for the contribution of «/ in the 
computation of C. In the case that R=\, then l-K", | = 1 and the contribution of each a/ is 

equal and therefore, to compute C with accuracy of three significant digits requires at 
least the five leading «/ as seen in Table 1 (the exact values); this can be seen by 
summing the exact values of Table 1, multiplied with the appropriate Ki (Ki =1 or -1). 
However, this is the extreme case where R is equal to the size of the original 
computational domain. In the usual range of values of R the required Na for the 
accurate computation of C is restricted to less than five «/ as it will be seen below. 



Table 3 
Dependence of Ki on R. 





K, 


K2 


^3 


K, 


^5 


R=0.9 


0.9486 


-0.8538 


0.7684 


-0.6915 


0.6224 


R=O.S 


0.8944 


-0.7155 


0.5724 


-0.4579 


0.3663 


R=0.7 


0.8366 


-0.5856 


0.4099 


-0.2869 


0.2008 



R=0.6 


0.7745 


-0.4647 


0.2788 


-0.1673 


0.1003 


R=0.5 


0.7071 


-0.3535 


0.1767 


-0.08838 


0.04419 


R=OA 


0.6324 


-0.2529 


0.1011 


-0.04047 


0.01619 


R=0.3 


0.5477 


-0.1643 


0.04929 


-0.01478 


0.004436 


R=0.2 


0.4472 


-0.08944 


0.01788 


-0.003577 


0.0007155 


R=0.l 


0.3162 


-0.03162 


0.003162 


-0.0003162 


0.00003162 



The Ki values for different values of R (Table 3) are derived from the initially known 
form of the solution and can be used as a 'loose' criterion for the number of «/ needed 
for the good approximation of C with respect to R. For example, for 7? = 0.3 a choice 
of Na with relatively good accuracy would be Na=2, given that K3 = 0.04929 
compared to Ki = 0.5477. The above criterion is only indicative since the value of «/ 
that multiplies Ki is unknown. However, it should also be taken into account that the 
absolute values of «/ decrease when / increases (for the Motz problem). For example, 
even though Ki = 0.94 and ^5 = 0.62 for 7? = 0.9 are very close, the contribution in C 
of the first term in (23) is -380 while the contribution of the fifth term is -0.89. In 
realistic applications R should only be a small portion of the size of the original 
computational domain. However, for extensively analyzing the proposed method, R is 
increased up to 0.9. 

The exact dependence of C on i? for the Motz problem and with different values of Na 
is shown in Fig. 5. The vertical axis corresponds to the relative error of the computed 
capacitance with respect to the exact for the given R, defined as 

£ = |Q-C|/Qxl00% (25) 

where Cex is computed by (23) with the ten exact leading singular coefficients of the 
Motz problem (cf. Table 2). All of the curves that correspond to different values of Na 
follow the same trend, declining at small values of R then reaching a minimum and 
start climbing again, some of which reaching error values that exceed 10% — for Na 
=1 we don't observe the above behavior, or at least for the range of R that was 
examined (0.01-0.9). Examining the graph of Fig. 5 we can extract the optimal values 
of TVct for intervals of i? that minimize E. For increasing Na with increments of 1 from 
Na=l to Na=5 the optimal intervals of 7? in the same order are: 7?=0-0.05, /?=0.05-0.15, 

7?=0. 15-0.25, 7?=0.25-0.4, R=OA-. 

10' 

10' 
10° 

10"' 

10-^ 
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R 

Fig. 5. Relative error, E, vs. radius of Qi, R; N=500, M=100. 




The proposed method is tested against the standalone BEM and the Galerkin/flnite 
element method (GFEM) [18, 19] in terais of the computed capacitance of the 
Dirichlet boundary, u of Fig. 4. On Fe, Cp^ is computed from (23) and on F5 

since the derivative of the solution is approximated via (4) (linear basis functions), a 
standard trapezoidal integration is adequate for the computation of Cp . The total 

capacitance is then C = Cj- + Cp . For the standard BEM and GFEM only the 

trapezoidal rule is applied on the whole lower left boundary. The computed C is 
compared with the exact C on Fj uFg computed from (23) with iV„=10. The relative 

error is then computed from (25), with Ce;c=340.30. For the various methods to level 
with each other, we chose a uniform discretization on the boundaries for the BEM and 
the proposed method, and uniform structured mesh for the GFEM with equally 
discretized boundaries; e.g. a 40x40 uniform GFEM mesh corresponds to 40 elements 
per side of Q that corresponds to N = 160 for the BEM. The current GFEM employs 
bilinear basis functions and thus, the 40x40 mesh gives 1681 degrees of freedom 
(DOF), while the BEM gives 160, yet for both of them the number of elements on F5 
of Fig. 4, denoted 7Vp_ is equal to 20; Te does not exist when the BEM or GFEM is 

applied. The difference between GFEM and BEM in terms of C on T5 is relatively 
small, with E following the same pattern with respect to A^p (Fig. 6). The GFEM 

manages to attain a value of E -3.2% for A'^p^ = 300 that corresponds to DOF=361201 

while the BEM attains the same value for Nr = 400 and DOF=3200 that also 

corresponds to A^=3200; however, the increased DOF of the GFEM is compensated by 
the sparsity of its matrices, while their counterparts of the BEM are dense. For the 
comparison of BEM and GFEM with the proposed hybrid method, three different 
values of R are used that are accompanied by the optimal Na that minimize E (cf. Fig. 
5). Moreover, as in the BEM method the hybrid method uses, 
TV = 4A^P^^P +M = SA'^P^ +M and for further simplicity, M is constrained as, 

M = 27Vp^ and thus, DOF=A^ + M + A^„ =12A^r^ +A^„. For all values of R the 

proposed method even for small TVp achieves small values of E, -1-2%. This is in 

contrast with the results from Fig. 5 (which are even smaller), due to C in Fig. 5 being 
computed only on Fg while in Fig. 6, C is computed on u Fg . 

We can safely assume that the discrepancies between Cp^ and Cp^^p^ are related to the 

integration on Fg, i.e. the boundary treated by the BEM. Moreover, the decrease ofE 
of the proposed method with respect to the increasing R is due to the decrease of F5 
and thus, the contribution of F5 on the total capacitance. The miscalculation of the 
integral on F5 is largely affected by the oscillations of the solution qj close to the ends 
of F5 (Fig. 7), which is a common behavior of the BEM close to comers of domains or 
at the points where the boundary conditions change. To further illustrate the effect of 
the oscillations of q on Cp^ and thus, on the total C, we seek to eliminate them by 

simply excluding the solution near the ends of F5 and fitting a curve based on a 
nonlinear regression model applied on the remaining values of qj (Fig. 7). The integral 
of the fitted curve is then the corrected Cp . This is done at the post-processing stage 

and does not belong to the core of the proposed method, but is done only to justify the 
above remark. However, there are techniques for the BEM, e.g. discontinuous 
elements at the corners, which provide more accurate solution near the comers, but 
are not incorporated in the present work. 
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Fig. 6. Relative error, E, vs. N ^-^ ; comparison between the GFEM, BEM and hybrid integral method. 

The relative error, E, based on corrected capacitance, Cp^ , which in turn is computed 

from the integral of the fitted curve of Fig. 7, is presented in Fig. 8. There is a clear 
reduction of E to levels that are seen in Fig. 5. This should provide sufficient proof 
that E computed on Fj u Fg is relatively large due to the oscillations of the solution, 
qj, on Fs near its corners. 




Fig. 7. Solution qj of the proposed method on Fs (solid line) and fitted curve (dots); =100. 
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Fig. 8. Relative error, E, vs. A'^p^ ; solid lines are identical to the ones in Fig. 6 (hybrid method) and 
dotted lines are derived from (25), implementing the corrected C that is computed by the fitting of on 

The proposed method is tested on a computational domain with two singularities. The 
problem is again the Motz problem even though it contains only one prominent 
singularity. However, every comer of the computational domain of the Motz problem 
can potentially impose a singularity that varies in strength, but comparatively to the 
singularity at (0, 0) is negligible. Despite this, we treat the area at the upper left comer 
as we treated the area at (0, 0). The computational domain is decomposed as in Fig. 9. 
Homogeneous Neumaim-Neumann conditions are imposed on the straight segments 
of the boundary of Q.2 and the asymptotic expansion of the solution is given by (1) 
with = nil and /// = 2(/ -1) [1]. 

The subdomain does not contain a singularity in a tme sense because the 
derivative of the solution in the radial direction does not reach infinity when r reaches 
0; this can be seen by differentiating (1) with r and using //^ = 2{k - 1) . However, the 

goal of this experiment is to analyze the behavior of the proposed method with the 
original computational domain decomposed into three subdomains as in Fig. 2. The 
parameters of this experiment that are fixed are A^= 400 and M= 100, while R\ and Ri 
are varied from 0.05 to 0.5 and from 0.1 to 0.4, respectively. The above constitutes a 
set of experiments for a given N ^ ^^^ . In Fig. 10 are presented three sets of 

experiments for N^^^^ = N^^^^ = 2,3,5 with their corresponding curves clustered 

around the dashed curve, which in turn corresponds to the same problem setup (i.e. 
input parameters) but applied on the computational domain of Fig. 4 (one singularity 
only); it should be noted that for these set of experiments we employed constant 
elements instead of linear. Each cluster of curves corresponds to a different A^^ and 

each curve from each cluster corresponds to a different R2. The deviation of the solid 
lines from the dotted refiects the effect of the presence of ^2 on the capacitance error 
of Fe, defined by (25). As it can be seen from Fig. 10, the effect of Q.2 for A^^ = 2 is 

minimal throughout the whole range of i?i and R2. For A^„ n, = 3 the effect is more 

prominent only at the range of R\ ~0. 15-0.25, while for A^„ q_ =5 the effect is 

prominent at the range of Ri -0.25-0.5. It is under investigation whether the 
discrepancies can be attributed to certain aspects of the proposed method, such as the 
proximity of the two singularities or the ratio of the radii of the two subdomains. 




when the involved aspects are quite numerous. These discrepancies may come from 
outside the proposed method; e.g. the accuracy of the system solver can be 
compromised by the condition number of the A matrix when A^^ increases. 
From the resuhs of Fig. 10 we can deduce, however, that small values of Ri and R2 
renders the proposed method less sensitive to the other input parameter (No). This 
behavior of the proposed method is desirable in computational practice since there is 
no need to seek optimal values for several input parameters. Thus, the choice of small 
Ri and R2, combined with the fact that a small R requires only small Na (cf. Fig. 5), 
gives a general guideline for the preferable choice of the input parameters. 
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Fig. 9. Domain decomposition of the Motz problem with three subdomains. 
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Fig. 10. Relative error, E, vs. Ri; N=500, M=100, 7?2=0.1-0.4. 



4 Conclusions 



The hybrid BEM/SFBIM treats elliptic boundary value problems with singularities 
efficiently, providing adequate accuracy in the results while maintaining low 
computational cost in terms of domain discretization, degrees of freedom, etc. That is 
in contrast with the standalone BEM or GFEM, where accuracy comes at the expense 
of the computational efficiency. The proposed method can be seen as an augmented 
BEM with singular functions and it is seamlessly embedded in an existing BEM code. 
The method is best suited for applications that require the computation of the leading 
singular coefficients (generalized stress intensity factors) or the accurate computation 
of the derivative of the solution near singularities. The results are evaluated in terms 
of the capacitance of the Dirichlet boundary adjacent to the singularity and are 
compared with those of the BEM and GFEM. 
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